home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / rs_test.pro < prev    next >
Text File  |  1997-07-08  |  5KB  |  151 lines
  1. ;$Id: rs_test.pro,v 1.5 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       RS_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function tests the hypothesis that two sample popultions,
  11. ;       {X[i], Y[i]}, have the same mean of distribution against the
  12. ;       hypothesis that they differ. The result is a two-element vector
  13. ;       containing the nearly-normal test statistic Z and the one-tailed
  14. ;       probability of obtaining a value of Z or greater. This type of 
  15. ;       test is often refered to as the Wilcoxon Rank-Sum Test or Mann-
  16. ;       Whitney U-Test.
  17. ;
  18. ; CATEGORY:
  19. ;       Statistics.
  20. ;
  21. ; CALLING SEQUENCE:
  22. ;       Result = RS_test(X, Y)
  23. ;
  24. ; INPUTS:
  25. ;       X:    An n-element vector of type integer, float or double.
  26. ;
  27. ;       Y:    An m-element vector of type integer, float or double.
  28. ;
  29. ; KEYWORD PARAMETERS:
  30. ;      UX:    Use this keyword to specify a named variable which returns
  31. ;             the Mann-Whitney statistic for X.
  32. ;
  33. ;      UY:    Use this keyword to specify a named variable which returns
  34. ;             the Mann-Whitney statistic for Y.
  35. ;
  36. ; EXAMPLE:
  37. ;       Define the vectors of sample data.
  38. ;         x = [-14,   3,   1, -16, -21,   7,  -7, -13, -22, -17, -14, -8, $
  39. ;                7, -18, -13,  -9, -22, -25, -24, -18, -13, -13, -18, -5]
  40. ;         y = [-18, -9, -16, -14,  -3,  -9, -16, 10, -11, -3, -13, $
  41. ;              -21, -2, -11, -16, -12, -13,  -6, -9,  -7, -11, -9]
  42. ;       Test the hypothesis that two sample popultions, {X[i], Y[i]},
  43. ;       have the same mean of distribution against the hypothesis that they
  44. ;       differ at the 0.05 significance level.
  45. ;         result = rs_test(x, y, ux = ux, uy = uy)
  46. ;       The result should be the 2-element vector:
  47. ;         [1.45134, 0.0733429]
  48. ;       The keyword parameters should be returned as:
  49. ;         ux = 330.000, uy = 198.000
  50. ;       The computed probability (0.0733429) is greater than the 0.05
  51. ;       significance level and therefore we do not reject the hypothesis
  52. ;       that X and Y have the same mean of distribution. 
  53. ;
  54. ; PROCEDURE:
  55. ;       RS_TEST computes the nonparametric Rank Sum Test for populations of
  56. ;       equal or unequal size. The populations X[i] and Y[i] are combined
  57. ;       and individual elements are ranked based on magnitude. Elements of
  58. ;       identical magnitude are ranked using a rank equal to the mean of the
  59. ;       ranks that would otherwise be assigned. The Mann-Whitney statistics
  60. ;       (Ux and Uy) are computed and used to determine the nearly-normal test
  61. ;       statistic (Z) and the one-tailed probability of obtaining a value of 
  62. ;       Z or greater. The hypothesis that two sample populations have the same
  63. ;       mean of distribution is rejected if Ux and Uy differ with statistical 
  64. ;       significance. If either population contains 10 or fewer samples, the
  65. ;       test statistic (Z) and the one-tailed probability of obtaining a value
  66. ;       of Z or greater are returned as zero. In this case, consult published
  67. ;       tables such as the ones available in the REFERENCE given below. 
  68. ;
  69. ; REFERENCE:
  70. ;       PROBABILITY and STATISTICS for ENGINEERS and SCIENTISTS (3rd edition)
  71. ;       Ronald E. Walpole & Raymond H. Myers
  72. ;       ISBN 0-02-424170-9
  73. ;
  74. ; MODIFICATION HISTORY:
  75. ;       Written by:  GGS, RSI, August 1994
  76. ;-
  77.  
  78. pro crank, w, s
  79.   ;Replace elements of the sorted array "w" by their rank.
  80.   ;Identical observations ("ties") are ranked according to their means.
  81.   ;s = f^3 - f (f is the number of elements in identical observations.)
  82.   n = n_elements(w)
  83.   w = [0.0, w]  ;operate on elements w[1], ... , w[n] of the shifted
  84.                 ;n+1 element float array (w).
  85.   s = 0.0
  86.   j = 1
  87.   while(j lt n) do begin
  88.     if(w[j+1] ne w[j]) then begin
  89.       w[j] = j
  90.       j = j+1
  91.     endif else begin
  92.       for jt = j+1, n do $
  93.         if (w[jt] ne w[j]) then goto, case2
  94.       jt = n + 1
  95.       case2:
  96.       rank = 0.5 * (j + jt - 1)
  97.       for ji = j, jt-1 do $
  98.         w[ji] = rank
  99.       t = jt - j
  100.       s = s + t^3 - t
  101.       j = jt
  102.     endelse
  103.   endwhile
  104.   if(j eq n) then w[n] = n
  105.   w = w[1:*]
  106. end
  107.  
  108. function rs_test, x, y, ux = ux, uy = uy
  109.  
  110.   on_error, 2
  111.  
  112.   nx = n_elements(x)
  113.   ny = n_elements(y)
  114.  
  115.   ;An alternate method of error handling if Ux and Uy are not required.
  116.   ;if min([nx, ny]) le 10 then message, $
  117.   ;    'Consult statistical tables for samples of 10 or fewer elements.'
  118.  
  119.   ;Number of "ties" (identical data).
  120.   ties = where(x - y eq 0, zdiff)
  121.  
  122.   if nx eq ny and zdiff eq nx then message, $
  123.     'x and y contain identical data.'
  124.  
  125.   ;Sort and rank the combined x and y vectors.
  126.   xy = [x, y]            ;Combine vectors.
  127.   is = sort(xy)          ;Sort xy into ascending order.
  128.   xy = xy[is]
  129.   crank, xy              ;Rank sorted xy vector.
  130.   xy = xy[sort(is)]
  131.   
  132.   ;Rank-sum totals.
  133.   wx = total(xy[0:nx-1])
  134.   wy = total(xy[nx:*])
  135.  
  136.   ;Compute the Mann-Whitney statistics.
  137.   ux = nx * ny + (nx * (nx + 1.0) / 2.0) - wx
  138.   uy = nx * ny + (ny * (ny + 1.0) / 2.0) - wy
  139.  
  140.   ;Compute the Z statistic with respect to ux.
  141.   z = (ux - (nx*ny/2.0)) / sqrt(nx*ny*(nx+ny+1.0)/12.0)
  142.  
  143.   ;Probability of obtaining z or something greater.
  144.   prob = 1.0 - gauss_pdf(abs(z)) 
  145.  
  146.   ;If either sample is <= 10, consult published statistical tables.
  147.   if min([nx, ny]) le 10 then return, [0, 0] $
  148.   else return, [z, prob] 
  149.  
  150. end
  151.